Astronomy & Astrophysics manuscript no. art 


©ESO 2012 


September 28, 2012 





Periodic morphological changes in the 
radio structure of the gamma-ray binary LS 5039 

J. MoldorfD M. Ribfl and J.M. ParedelO 

'Departament d'Astronomia i Meteorologia, Institut de Ciencies del Cosmos (ICC), Universitat de Barcelona (IEEC-UB), Marti i 
Franques 1, 08028 Barcelona, Spain 

e-mail: jmoldon@am.ub.es; mribo@am.ub.es; jmparedes@ub.edu 

Received / Accepted 

ABSTRACT 

Context. Gamma-ray binaries allow us to study physical processes such as particle acceleration up to TeV energies and very-high- 
energy gamma-ray emission and absorption with changing geometrical configurations on a periodic basis. These sources produce 
outflows of radio-emitting particles whose structure can be imaged with very long baseline interferometry (VLBI). LS 5039 is a 
gamma-ray binary that has shown variable VLBI structures in the past. 

Aims. We aim to characterise the radio morphological changes of LS 5039 and discriminate if they are either repeatable or erratic. 
Methods. We observed LS 5039 with the VLBA at 5 GHz during five consecutive days to cover the 3.9-day orbit and an extra day to 
disentangle between orbital or secular variability. We also compiled the available high-resolution radio observations of the source to 
study its morphological variability at different orbital phases. We used a simple model to interpret the obtained images. 
Results. The new observations show that the morphology of LS 5039 up to projected distances of 10 milliarcseconds changes in 
24 h. The observed radio morphological changes display a periodic orbital modulation. Multifrequency and multiepoch VLBI obser- 
vations confirm that the morphological periodicity is stable on timescales of years. Using a simple model we show that the observed 
behaviour is compatible with the presence of a young non-accreting pulsar with an outflow behind it. The morphology is reproduced 
for inclinations of the orbit of 60-75°. For masses of the companion star in the range 20-50 M Q , this range of inclinations implies a 
mass of the compact object of 1.3-2.7 M Q . 

Conclusions. The periodic orbital modulation of the radio morphology of LS 5039, stable over several years, suggests that all gamma- 
ray binaries are expected to show a similar behaviour. The changes in the radio structure of LS 5039 are compatible with the presence 
of a young non-accreting neutron star, which suggests that the known gamma-ray binaries contain young pulsars. 

Key words, stars: individual: LS 5039 - radio continuum: stars - binaries: close - gamma rays: stars - X-rays: binaries - radiation 
mechanism: non-thermal 



1. Introduction 

Gamma-ray binaries are extreme systems that produce non- 
thermal emission from radio to very-high-energy (>TeV) 
gamma rays, with the energy output in the spectral energy 
distribution (S EP) dominated by the MeV-GeV photons (see 
1 lParedesl201 ll for a review). Their broadband emission is usually 
modulated by the orbital cycle of the system, which suggests that 
the physical conditions are also periodic and reproducible. The 
wide range of different orbital periods and eccentricities of th e 
known gamma-ray binaries (see Table 4 in lCasares et alJ2012bl) . 
provides a diversity of different ambient conditions in which the 
physical processes take place. The diversity of systems, together 
with the reproducibility of the conditions within each system, 
makes gamma-ray binaries excellent physical laboratories in 
which high energy particle acceleration, diffusion, absorption, 
and radiation mechanisms can be explored. Nevertheless, 
the number of known gamma-ray binaries is still very lim- 
ited. Five binary systems have been classified as gamma-ray 
binaries and present orbital modulation of the GeV and/or 
TeV gamma-ray emission: PSR B1259-63 ([A haronian et alj 
l2009h . LS 5039 (lAharonian etail l2006t lAbdo et all |2009b). 
LS I +61 303 (lAlbert et alJ 120091: lAbdo et alJ Ep )09a), 
HESS J0632+057 ( Aharoni anet alJ 120071: iMaier et all l2oTTb 
and 1FGL J1018.6-5856 dFermi LAT Collaboration et alJ 



120121) . The latter has not been clearly detected at TeV energies 
yet. On the other hand, Cygnus X-3 presents active periods 
of GeV emission in which o r bital pe riodicity is d etecte d 
dFermi LAT Collaboration etaf] l2009t iTavani et alJ [2009). 
Cygnus X-l showed evidence of gamma-ra y emission at TeV 
energies on a single night of observations dAlbert et al. 2007) 
. However, in these two cases the peak of the SED is at keV 
energies and therefore are classified as X-ray binaries, the 
high energy activity of which is clearly powered by accre- 
tion. Another phenomenological difference between X-ray 
and gamma-ray binaries is that variations in the non-thermal 
broadband emission in X-ray binaries is usually linked to X-ray 
state changes, whereas in gamma-ray binaries it is coupled to 
the orbital period. 

LS 5039 is a binary system composed by the homonym 
bright star LS 5039, of spectral type ON6.5 V((f)), and a com- 
pact object dCasares et al.ll2005l) . The orbital period of the sys- 
tem is 3.9 days and the eccentricity of th e orbit is e ~ 0.3 5 
dCasares et alJ 120051: lAragona et all l2009t ISartv et all l201ll) . 
Given the mass function of the binary system the compact object 
would be a black hole for low inclinations of the or bit, and a neu- 
tron s tar for high inclinations (s ee the discussion in | Casares et al.l 
2005 and recent constraints in ICasares et alj|2012ah . No short- 
period pulsations were found that could demonstrate th e pres- 
ence of a pulsar either in radio (Mc Swain et alJ [201 ll) or X- 
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rays (iMartocchia et all 120051: iRea et alJ 120111) . The distance 
to the system has recently been updated to 2.9 + 0.8 kpc 
dCasares et alj2012ah . LS 5039 displays non-ther mal, persistent, 
and periodic gamma-ray emi s sion up to 4 TeV dParedes et alJ 
2000; lAharonian et all 120051: lAbdo etal.1 l2009bh . The X-ray 
and g am ma-ray emission show a p eriodic modulat i on of 
3.9 days dBosch-Ramon et alJ l2005t lAharonian et alJ l2006t 
iKishishita et al.ll2009HAbdo et al.ll2009bl) . The total radio emis- 
sion of LS 5039 is variable, altho ugh no periodicity or strong 
outbursts ha ve been detected so far (M artf et alJl 998: Ribo et all 
[19991120021) . The radio emission above 1 GHz is non-thermal, 
with spectral index of -0.46 (Marti et al. 1998 ), and inverted 
at low er frequ encies dGodambe et all 2008: Bhattacharvva et all 
120121 but see iPandev et al.ll2.Q07h . The synchrotron radio emis- 
sion appears extended when observed with VLBI on scales of 
milliarcseconds (mas). The source shows a main core and ex- 
tended bipolar emission that has been observed in directions 
with position angles (RA.) between 120 and 150° and projected 
angular distance between 1 and 180 mas ( 3-500 AU) from the 
core dParedes et alj|2000l 120021: iRibo et al.ll2008l) . In particular, 
Ribo et alJ d2008l) measured a change in the source morphology 
at scales of 2-6 m as from VLBI image s obtained 5 days apart. 
Finally, for X-rays, iDurant et alJ d201 lb discovered an extended 
component in X-rays up to 2' from LS 5039. 

Other gamma-ray binaries also present extended and vari- 
able radio emission at mas scales. PSR B1259-63 is the only 
gamma-ray binary in which radio pulsations have been de- 
tected up to now. During the periastron passage of this eccen- 
tric system, occurring every 3.4 yr, the close interaction be- 
tween the pulsar and the massive Be star and/o r Be disc pro- 
duces b roadband emission during a few months. Moldon et al. 
d2011al) discovered transient extended emission during and af- 
ter the periastron passage with a total projected extent of ~ 
50 mas (120 AU). The peak of the radio emission was de- 
tected outside the binary system. VLBI observations at differ- 
ent orbital phases of the 26.5-d period binary LS I +61 303, 
which also hosts a Be star, were conducted by IDhawan et alJ 
d2006h . showing an orbital variability that these authors inter- 
preted as the signature of a cometary tail produced in the collid- 
ing winds scenario. Another interpretation of the data suggests 
that the changes a re compatible with a precessing microblazar 
dMassi etalj|2012l) . The system HESS J0632+ 057, also hosting 
a Be star, has an orbital period of 321 days dBongiorno et alJ 
and faint radio extended emission at mas scales was de- 
tected ~ 130 days after the periastron passage (Moldon et al. 
I2011bl) . For 1FGL J1018 . 6-5856, with a period of 16 d and als o 
hosting an 06V((f)) star dFermi LAT Collaboration et alj|2012l) . 
no VLBI observations have been reported so far. In addition, 
some of these systems also present extended X-r ay emission at 
larger scales of ~ 1', as hinted for LS I +61 303 (|Paredes et ail 
2007) and shown for PSR B 1259-63 dPavlov et al.ll201 ll) . 

Several theoretical models were developed to explain the 
multiwavelength behaviour of LS 5039. The very-high-energy 
gamma-ray emission can be interpreted as the result of inverse 
Compton upscattering of stellar UV photons by relativistic elec- 
trons. The acceleration of electrons can be explained by two ex- 
clusive scenarios: acceleration in the jet of a microquasar pow- 
ered by accretio n dParedes et al. 2006; Bosc h-Ramon et aD2006l 
and the review in Bosch-Ramon & Khang ulvanl2 009). or shocks 
between the relativistic wind of a y oung non-accreting pulsar 
and the wind of the ste ll ar companion dMaraschi & Treves 1981; 
iTavani & Aroiisl 11997b iDubusf 120061: Ixhangulvan et alJ 120071) 
A simple and shockless microquasar scenario is disfavoured 
by previous VLBI radio observations of LS 5039 (Rib o et all 



2008). Also, no signs of the presenc e of an accretion disc 
have been detected so far (although see iBarkov & K hangulvan 
20121 and lOkazaki et alJ 120081 for alternative explanations). 



Sartv et all (1201 ll) used the stability of the optical photome- 
try of LS 5039 to constrain the orbit inclination to be be- 
low 60° and the mass of the compact object above 1.8 M Q . 
Models based on the young n on-accreting pul s ar sce nario, 
first proposed for LS 5039 in IMartocchia et alJ d2005l) . had 
provided descriptions of the HE/VHE light curves and the 
spectral evo lution of the source as a function o f the orbital 
phase (see ISierpow ska-Bar tosik & Torres! 120071 : Dubus etahl 
2008; iKhangulvan et alJ 120081: iTakahashi et alJ 120091) . Apart 
from the acceleration of particles in the main shock be- 
tween winds, additional gamma-r ay emission c an be produced 
by the unshocked pu lsar wind (|Cerutti et alJ 20081) and by 
secon dary cascading dBosch-Ramon et alT 20081 : ICerutti et al.l 
2010). The X-ray light curve and spectrum show orbital vari- 
ability, although there are no signatures of varia ble X-ray 
absorp t ion or X-ray occ ultatio ns, as discussed in iReig et all 
d2003[).lMartocchia et al1d2005h.lBosch-Ramon et al.ld2007l) and 
ISzostek & Dubusl d201 ll) . IZabalza et alJ d201 lb constrained the 
spin-down luminosity of the putative pulsar to be ~ 3-6 x 
10 36 ergs _1 . 

Very recently, detailed hydrodynamical simulation have 
been obtained to understand the wind shocks and the flow 
dynamics in gamma-ray binaries. Two-dimensional hydro- 
dynamic simulations were performed to study the wind- 
wind collision intera ction at scales of the binar y system 
for PSR B125 9-63 dBogoyalov et all 120081 [2012h and for 
LS I +61 303 dRomero et alJl2007l). Also fo r PSR B 1259-63, 
lOkazaki et alJ d201 ll) and Takata et alJ d2012|) presented 3D sim- 
ulations of the tidal and wind interactions. lLamberts et al.l d2012l) 
considered the general case of colliding wind binaries and de- 
scribed the outflow structure at larger distances. Finally, the 
flow structure for a case s imilar to LS 5039 was described in 
Bosch-Ramon et a ll d2012h for scales up to 100 times the orbit 
size. 

From a theoretical point of view, a general discussion on 
th e properties of the radio em ission of LS 5039 can be found 
in Bosch-Ram on et al.l d2008l) . However, only a few predi ctions 
on the expected radio morphology have been discussed. Dubus 
(2006) presented emission maps at different orbital phases for 
LS 5039 and predicted orbital morphological changes at mas 
scales, as well as displacements of the peak of t he emission of 
a few mas. iBosch-Ramon & Khangulvanl d201 ll) obtained syn- 
thetic maps at radio wavelengths from a model based on sec- 
ondary particles created from gamma-ray absorption, which may 
account for a significant fraction of the total radio flux density 
of the source and also predicts extended radio emission at mas 
scales. 

In this paper we present the results from a dedicated VLBA 
campaign that covers a full orbital cycle and provides astromet- 
ric and morphological information during five consecutive days. 
The observations and data reduction are presented in Sects. |2]and 
[3] including a discussion on the calibration caveats. In Sect.|4]we 
present the obtained astrometry and morphology of the source at 
different orbital phases and analyse the components character- 
ising the extended emission of the source. In Sect. [5] we show 
a compilation of VLBI data on LS 5039 and we describe and 
compare the morphology at different orbital phases for observa- 
tions at frequencies between 1.7 and 8.5 GHz taken from 1999 
to 2009. In Sect. [6] we present a model to check if the measured 
changes are compatible with a wind-wind collision scenario. The 
model allows us to estimate the orientation of the orbit on the 
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Table 1. Log of VLBA observations of project BR127, con- 
ducted at 5 GHz. The orbital phases have an uncertainty of 0.02. 



Run Date (Y-M-D) 



MJD 



Phase range 



A 


2007-07-05 


54286.15- 


-54286.40 


0.00-0.07 


B 


2007-07-06 


54287.15- 


-54287.40 


0.26-0.32 


C 


2007-07-07 


54288.15- 


-54288.40 


0.51-0.58 


D 


2007-07-08 


54289.15- 


-54289.40 


0.77-0.83 


E 


2007-07-09 


54290.15- 


-54290.40 


0.03-0.09 



sky, in particular its inclination, and consequently constrain the 
mass of the compact object. Finally, we present a summary of 
the results and the main conclusions of this work in Sect. [7] 



2. Observations 

We conducted VLBI observations on LS 5039 during five con- 
secutive days to cover an orbit of 3.9 d and an extra day to 
disentangle between orbital or secular variability. We observed 
at 5 GHz (6 cm wavelength) with the Very Long Baseline 
Array (VLBA) of the National Radio Astronomy Observatory 
(NRAO). The stations forming the array are Br, Fd, Hn, Kp, 
La, Mk, Nl, Ov, Pt, and Sc. The VLBA project code is BR127, 
and the five observations were conducted from July 5 to 9, 
2007 (MJD 54286 to 54290). The five runs, hereafter run A- 
E, spanned from 03:30 to 09:30 UTC and were scheduled in 
the same way to reduce differences between runs. The corre- 
sp onding orbital phases were computed using the ephemerides 
in lCasares etall d2012ah . T = HJD 2453 478.09(6) and P orb = 
3.90603(8), where the values in parentheses refer to the uncer- 
tainty in the last digit. The phase of the periastron passage is 0.0. 
The orbital phase at the centre of each observation was 0.03, 
0.29, 0.55, 0.80, and 0.06, respectively, with an uncertainty of 
0.02, or approximately 1.9 hours. A log of the observations is 
shown in Table Q] 

The data were obtained with single circular left-handed po- 
larisation, with eight sub-bands of 8 MHz, and were correlated 
with two bits per sample, obtaining a total aggregate bit-rate of 
256 Mbps. The data were processed at the VLBA hardware cor- 
relator in Socorro, which produced the final visibilities with an 
integration time of 2 seconds. 

We conducted the observations using the phase-referencing 
technique on the phase reference calibrator J 1825- 17 18, located 
at 2.5° from LS 5039 in a P.A. of -176.4° (see Fig. Q). We 
observed in 4.2-min cycles, spending 2.5 minutes on the tar- 
get source and 1.7 minutes on J1825-1718. The total scheduled 
time on the target source was 2.7 hours. As an astrometric check 
source we used J1818— 1 108, including one scan of 2.5 min- 
utes every 30 minutes. Finally, we included three 5-minute scans 
on the fringe finder J1733-1304. The reference position for the 
observations is the correlation position of the phase reference 
source J1825-1718 a J2 ooo.o = 18 h 25 m 36 s .53228 ± 0:00009 (or 
+ 1.3 mas) a nda.i?nno.o = -17°18 , 4 9 / /848 + 07002 (or +2 mas), 
taken from iFomalont et alj (2003) via the SCHED database 
(NASA catalogue 2005f_astro). LS 5039 is close to the galac- 
tic plane, at a galactic latitude of -1.29°. The phase calibra- 
tor and the astrometric check source suffer from galactic scat- 
te r broadening (a g eneral discussion on this effect can be found 
in lFev et alJll99lh . The phase calibrator J1825-1718 has a to- 
tal flux density of 350 mJy at 5 GHz, but the amplitude of the 
visibilities decreases with the uv distance and the phase cali- 



bration fails for baselines longer than ~90-100 MA. The com- 
pact core has a size of 3 mas. The resolution of the phase- 
referenced images is therefore limited by the scatter broadening 
of the phase calibrator. On the other hand, the astrometric check 
source J1818— 1 108 is not compact. The self-calibrated images 
show a main core of 680 mJy and a secondary blob of 180 mJy 
located at 50 mas westwards from the core. Both components 
have a size of ~8 mas. The visibility amplitudes quickly drop 
with the uv distance and there is no signal beyond ~30 MA. The 
distance between J 1 8 1 8— 1 108 and the phase calibrator is 6.4° 
in RA. of 16° (see Fig. [TJ, which makes difficult to transfer the 
phase solutions for this low elevation source. 



-10° 



-12° 



-14° 



-16° - 



-20° 



J1818-1108 




J1825-171 



Direct phase referencing 
Inverse phase referencing 



18 h 36 m 18 h 24 m 18 

Right Ascension 



12 n 



Fig. 1. Distribution on the sky of LS 5039, the phase calibra- 
tor J1825-1718, and the astrometric check source J 1 8 18— 1 108. 
The dashed line indicates the Galactic Plane. The arrows indicate 
the two phase calibrations obtained: direct, using J 1825- 17 18 as 
phase reference; inverse, using LS 5039 as phase reference. 



3. Data reduction 

Data redu ction was princip ally performed in AIPSQ- The Difmap 
package (Shepherd| 1 1997b was used for imaging and self- 
calibration. The raw visibilities were loaded into AIPS, and a 
priori flagging on telescope off-source times because of antenna 
slewing was applied. We searched for data with instrumental 
problems and flagged them, as well as all visibilities with an- 
tenna elevation below 5°. Initially, we updated the geometrical 
model of the correlator using the Earth orientation parameters 
(EOPs). However, the correction implied a change on the final 
measured positions one order of magnitude below our final un- 
certainties. Considering the small effect and the problems with 
this correction in the past (wrong parameters during the corre- 
lations between 2003 and 20050 and a bug in the AIPS task 
CLCOR from 2009 to 2011), we did not include them in the 



1 The NRAO A stronomical 
http://www.aips.nrao.edu/ 



Image Processing System. 



http://www.vlba.nrao.edu/memos/test/test69memo/eop_problem.html 
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final data processing. For these observations, the EOPs correc- 
tion produced a ~ 10° phase offset in the visibilities, which was 
in any case removed by the phase calibration (see below). 



3.1. Ionospheric correction 

The different ionospheric conditions above VLBI antennas are 
one of the main contributions to systematic errors for astrometry 
at low frequencies. The unmodelled phase and delay contribu- 
tion of the ionosphere modifies randomly the observed phases 
on timescales of minutes, producing two effects. On one hand, it 
spreads the signal and broadens the source, producing an over- 
all decrease of the detected signal-to-noise ratio (S/N). On the 
other hand, it introduces an unknown position offset. This cor- 
rection is specially important in this case because the sources 
have low declination, which implies low elevations during the 
observations, and because the sources are separated mainly in 
the North-South direction (see Fig. [TJ, where the atmospheric 
conditions change more significantly. 

We used ionospheric total electron content (TEC) models 
based on GPS data obtained from the CDDIS data archival to 
correct the phase variations caused by the ionosphere. These 
maps provide one TEC measures every two hours on a grid of 
5° x 2.5° in terrestrial longitude and latitude, respectively. This 
is a coarse correction because the ionosphere changes occur on 
a shorter timescale, although it accounts for significant astro- 
metric offsets. The ionospheric models are produced by differ- 
ent institutes, which estimate the TEC above different locations 
and provide electron distribution maps (IONEX files) that can 
be loaded into AIPS to apply the phase corrections. We com- 
pared the effects of the phase correction on the visibilities us- 
ing the models provided by four institutes: the Jet Propulsion 
Laboratory (JPL), the Center for Orbit Determination in Europe 
(CODE), the ESOC Ionosphere Monitoring Facility (ESA), and 
the Universitat Politecnica de Catalunya (UPC). We tested the 
models by comparing the S/N and the stability of the position of 
the peak of the emission of LS 5039 in different short intervals of 
1 and 3 hours along the observations. We did not measure appre- 
ciable differences in the peak flux density when using the differ- 
ent models. The ionospheric models provide an average position 
offset of 0.02 and -1.7 mas, in right ascension and declination, 
respectively. The models are very similar, and the dispersion of 
these corrections are 0.11 and 0.2 mas, respectively. For each run 
we applied the model that provides offsets closer to the averaged 
value. The JPL model was used to correct the runs A and C, and 
the CODE model to correct runs B, D, and E. 

We checked the stability of the correction for the observed 
times by computing the variability of the TEC provided by the 
GPS maps. The most variable stations were Sc (at the end of the 
observation the Sun was already rising) and Mk (the observa- 
tions started before sunset). Comparing the mean variability of 
the TEC content for every run during the central hours of the 
project (from 4:00 to 8:00 UTC), we found that runs B, D, and E 
have a relative average variability of 18%, for run C it was 15%, 
and for run A 22%. This shows that during run A the ionosphere 
was, in general, more variable. This can in part justify the astro- 
metric problems encountered for this run (see Sect. 13.31 1. being 
the rest of the weather conditions similar for all runs. 



3 The Crustal Dynamics Data Information System 
http://cddis.nasa. gov/ 



3.2. Amplitude and phase calibration 

ACCOR was used to fix the amplitudes in the cross correlation 
spectrum from the VLBA correlator. The amplitude calibration 
was performed using the antenna gains and the system temper- 
atures measured at each station in real time during the observa- 
tions, using the APCAL procedure. The parallactic angle cor- 
rection was automatically applied with VLBAPANG. To correct 
the instrumental offsets and slopes between and within the dif- 
ferent bands we ran FRING on 30 seconds of a scan of the fringe 
finder 3C 345 (manual phase-cal). To correct the dependence of 
the visibility amplitudes with the frequency we used the auto- 
correlation values from 3C 345 to smooth the bandpass shape of 
the amplitudes. 

An initial fringe fitting was performed using the FRING rou- 
tine in AIPS, using a point-source model on the phase calibra- 
tor J 1825- 17 18. These data were exported to Difmap and self- 
calibrated. We fitted the visibilities for each run with a circular 
Gaussian. The mean flux density fitted for the core component 
was 345 mJy, with a standard deviation between runs of 2.3 mJy 
(0.67%), while the calibrator structure did not change signifi- 
cantly. Therefore, we combined the five data sets of the calibrator 
and performed several iterations of imaging and self-calibration 
to create a master calibrator image. We ran FRING again, but 
now we removed the source structure contribution to the phase 
calibration by using the combined image of J 1825- 17 18 as an 
input model. One solution for each scan for all IFs was searched 
within a delay and rate search windows of 80 nanoseconds and 
20 mHz, respectively. The minimum allowed S/N for solutions 
was 8, providing 95.8% of good solutions. The 4.2% of the data 
without solutions mostly correspond to long baselines (above 
90 MA). The solutions for every frequency band were corrected 
with MBDLY, which fits the phase delay between bands. The fi- 
nal phases were explored and flagged when necessary. Finally, 
we applied the calibration and flag tables, and we averaged all 
frequency channels within each IF to obtain three single-source 
files, one for LS 5039, one for the phase calibrator J 1825- 17 18 
and one for the astrometric check source J 1 8 1 8— 1 108. 

We also obtained inverse phase referencing based on 
LS 5039 because the source can be self-calibrated at 5 GHz. We 
followed an analogue procedure to fringe-fit the LS 5039 data, 
although we did not produce a combined data set. The S/N cut 
was set to 5 because of the faintness of the source. We flagged the 
unstable phases and transferred the solutions from LS 5039 to 
J1825— 17 18 and J1818-1 108, for which inverse astrometry was 
obtained. A sketch of the phase referencing schemes is shown in 
Fig.UJ 

3.3. Phase referencing imaging 

To reduce the systematic errors of the phase referencing pro- 
duced by data taken when the sources were at low elevations, 
we only used the central hours of the runs. We produced im- 
ages with time intervals of 30 and 60 minutes and inspected the 
image quality and the astrometric stability. We also produced 
images for the individual scans. The data obtained at the begin- 
ning and at the end of the observations are not consistent because 
we measured artificial displacements of the peak of the emission 
between ~0.5 and 3 mas, while the source became elongated in 
random directions. This behaviour can be explained by the iono- 
sphere affecting the data with low elevation. We note that the ef- 
fect is more important for run A, which suffered severe displace- 
ments of more than 6 mas in the time blocks for the first hour and 
the last two hours. After inspecting these systematic errors for 
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different time blocks we only used the data in the time intervals 
when the source was reasonably stable, from 04:40 to 07:40 UT 
for runs B, C, and D, from 04:30 to 07:00 UT for run A, and from 
06:00 to 08:30 UT for run E. Taking into account that the dec- 
lination of the phase-reference calibrator and the target source 
are -17.3° and -14.8°, respectively, that the observations were 
centred at the culmination of the sources, and that the observa- 
tions were conducted under normal weather conditions (without 
rain, strong winds or snow), the time ranges used correspond to 
elevations of most of the antennas above ~ 25°. Some antennas 
contributing to the longest baselines (Mk, Sc, Hn, and Br) were 
below this value during part of the observations, and the phase 
calibration was not succesful for them during these time inter- 
vals, altough this can be also caused by the intrinsic calibrator 
structure. 

The phase-referenced data were imaged using task IMAGR 
with a pixel size of 0.2 mas. We used a weighting scheme with 
robust as a compr o mise b etween angular resolution and sen- 
sitivity. iPradel etal] d2006l) showed that removing the shorter 
baselines of the VLBA can decrease the astrometric systematic 
errors by a 15%. For LS 5039, the source was detected with S/N 
above 10 when imaging the baselines with uv distances between 
15-20 and 100 MA. The range of minimum baselines was low- 
ered in those cases in which secondary lobes were more impor- 
tant or when the S/N of the peak emission was too low. The 
same strategy was used to image the astrometric check source 
J1818— 1 108. In this highly resolved source we only used base- 
lines with uv distances between 12 and 40 MA, 

For the inverse phase referencing, where phases from 
LS 5039 were transferred to J1825-1718 and J1818-1108, a 
uv range of 20 to 60 MA was used for J 1825- 17 18 and 12 to 
60 MA for the check source. The maximum baseline of the lat- 
ter is higher than in the direct phase referencing because of the 
smaller angular distance between J 1 8 1 8 - 1 1 08 and the phase ref- 
erence source in this case (see Fig. Q]). In both cases we used a 
pixel size of 0.5 mas. 

3.4. Self-calibration 

In order to perform the self-calibration of the LS 5039 data, we 
loaded the individual data files into Difmap. We averaged the 
data with a binning size of 50 seconds. The visibilities with 
long baselines were down-weighted, using a Gaussian taper at 
radius of 80 MA. We performed several iterations of cleaning 
and phase self-calibration. For each of them, phase solutions 
were found, first for the most compact part of the core by imag- 
ing the data with a uniform weighting scheme and then solv- 
ing for the more extended emission using a natural weighting 
cleaning. After each cycle, an amplitude self-calibration was ob- 
tained, each time with a shorter integration time, up to a mini- 
mum of 30 minutes. Finally, two hybrid images were produced 
for each run, one with natural weight (worse resolution and bet- 
ter sensitivity) and one with uniform weight (better resolution 
and worse sensitivity). The final images were exported to AIPS. 
Their rms noise and the total flux density were obtained by fit- 
ting the pixel flux density distribution of the image with IMEAN. 
The extended emission was characterised by fitting two or three 
Gaussian components to the images using the task JMFIT. 

We also fitted Gaussian components to the calibrated uv data 
sets of LS 5039 (modelfit). This is an iterative process, where 
the initial values for the fit have to be fixed manually, and the fi- 
nal solution can depend on the starting parameters. Therefore, 
we note that the model fitting results can be slightly subjec- 
tive, in particular for this case, where the source structure is not 



clearly defined by individual components. Despite this caveat, 
the modelfit provides morphological information directly from 
the uv data and therefore it is not limited by a synthesized beam 
because it is independent of the deconvolution procedure used. 
We obtained the modelfit of the uv data with the task UVFIT in 
AIPS. As a starting point, we used the solutions from JMFIT ob- 
tained from the naturally weighted images. Then we inspected 
several combinations of number of components and shape re- 
strictions until a robust solution was found. As a parallel check, 
we performed another modelfit in Difmap to verify the consis- 
tency of the results. As both approaches provide similar results, 
we only present here the solutions obtained with AIPS. 

4. Analysis and results 

4.1. Astrometry 

The average position of LS 5039 from the five runs of project 
BR127 is Qrj2oooo = 18 h 26 m 15!06003(3) (or +0.4 mas) and 
<5j2uoo.o = -14°50'54'.'3094(6) (or +0.6 mas). This posi- 
tion is measured with respect to the correlation position of 
J1825-1718, which is listed in Sect. [2] The absolute astrome- 
try in this proj ect was used to ob tain an accurate proper motion 
of LS 5039 dMoldon et al.ll2012l) . The astrometric uncertainties 
cannot be obtained directly from the position fit because of the 
unknown ionospheric effects, and therefore we used the astro- 
metric check source J18 18— 1 108 to estimate them. The standard 
deviation of the peak positions of the check source was 2.3 mas 
in right ascension and 2.4 mas in declination. These deviations 
were converted to uncertainties in the determination of the posi- 
tion of LS 5039 by multiplying them by two correction factors 
that account for the different observing conditions for LS 5039 
and J1818— 1 108, following th e general theoretical a strometric 
precision for an interferometer (Thom pson et al.l ll986) 



where the wavelength (A) and the maximum baseline (B) can 
be represented by the synthesized beam size. First, the uncer- 
tainty was scaled taking into account the different resolutions 
of the images. Basically, the data of J 1 8 1 8— 1 108 were imaged 
with poorer resolution, only using short baselines (see Sect. l3.3l >. 
This factor was obtained as a ratio of the synthesized beam size 
of the LS 5039 observations to the mean size of the J1818— 1 108 
beams. Second, the astrometric precision depends on the S/N of 
the detection. We scaled the uncertainty by the ratio of individ- 
ual S/N of LS 5039, which was around 15, to the mean S/N of 
J1818— 1 108, which was 7. These two corrections yield a scal- 
ing factor applied to the dispersion of J 1 8 1 8— 1 108 of ~ 0.19 and 
~ 0.25 in right ascension and declination, respectively. Finally, 
the mean uncertainties in the determination of the right ascen- 
sion and declination of LS 5039 are 0.4 and 0.6 mas, respec- 
tively. The same procedure was used for the inverted astrometry. 

Another approach to estimate the astrometric uncertainties 
is to use the expect ed theoretical syste matic uncertainties com- 
puted using Eq. 2 in lPradel etal] (2006). although some assump- 
tions considered in the general discussion in that paper are not 
met in these observations. The minimum theoretical uncertain- 
ties for LS 5039 using J 1825- 17 18 as a reference source and 
assuming mean tropospheric conditions are 0.17 and 0.56 mas 
in right ascension and declination, respectively, which are com- 
patible with the estimation quoted above. This indicates that the 
uncertainties in declination seem to be dominated by system- 
atic uncertainties. On the other hand, the minimum theoretical 
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uncertainties for J 1 8 1 8 — 1 108 are 0.43 and 1.42 mas in right as- 
cension and declination, respectively. The real dispersion mea- 
sured, 2.3 mas in right ascension and 2.4 mas in declination, is 
significantly larger possibly because the extended nature of the 
source, the poor uv coverage, and the low S/N. Given these re- 
sults, we have adopted the more conservative approach described 
above, although we note that the real uncertainties of the mea- 
surement are possibly a non-trivial combination of the system- 
atic and nominal uncertainties. 

In Table|2]we present the measured astrometry. We show the 
relative astrometry (with respect to the average position) for the 
check source J1818— 1 108 and the synthesized beam of the phase 
referenced images and the relative astrometry for LS 5039. We 
also show the inverted astrometry obtained using LS 5039 as 
a phase reference source. The measured relative offsets of the 
peak of the emission of LS 5039 using the two methods are plot- 
ted in Fig. |2 Both approaches provide similar relative astrome- 
try. The peak position of LS 5039 during run A is significantly 
displaced from the mean value. However, in run A, the check 
source shows a even larger displacement in the same direction 
(see first row in Tabled. This fact, combined with the suspicious 
behaviour described in Sect. 13.31 makes clear that the astrometry 
of run A suffers from uncorrected ionospheric effects, and there- 
fore we do not consider the peak displacement in run A as reli- 
able. The only significant displacement is found between run B 
and run C, which corresponds to displacements of 1 .0 + 0.5 mas 
and -2.0 + 0.7 mas in right ascension and declination, respec- 
tively. Therefore, the total displacement between phase 0.29 and 
phase 0.55 is 2.3 ± 0.6 mas in P.A. of 154°, which, considering 
the limitations of the data, cannot be considered a robust sign of 
peak displacement. 



Direct phase referencing 
Inverse phase referencing 



Orbit size 






B 



2 1 0-1-2 

R.A. offset [mas] 

Fig. 2. Relative offsets of the peak of the emission of LS 5039 
during five consecutive days with respect to the average position. 
For straightforward comparison, we plot the opposite value of 
the inverse astrometry. The run label is displayed close to each 
direct measure. The size of the binary system is indicated by 
the small orbit in the corner, which is plotted face-on and with 
an arbitrary orientation. The displacements and uncertainties are 
shown in Table [2] 



4.2. Morphology 

The self-calibration of the data improves the S/N of the phase- 
referenced images while preserving the main features of the ex- 
tended emission, partially avoiding the problems with the phase 
calibration caused by the ionosphere, at the expense of losing 
the astrometric information. In the first two rows of Fig. [3] we 
show the final self-calibrated images for each run, produced 
with natural and uniform weight, respectively. The five images 
with natural weight have synthesized beams with sizes about 
6.1 x 2.3 mas 2 and P.A. of -2° and total flux densities of 27.9, 
26.1, 24.4, 28.5, and 26.4 mJy, respectively, with a rms noise of 
the images of ~0.06 mJy beam 1 . On the other hand, the images 
with uniform weight have synthesized beams with sizes around 
4.1x1.3 mas 2 and P.A. of -4° . The total flux densities in the five 
images are 27.8, 26.7, 25.4, 28.5, and 27.0, with a rms noise of 
0.08 mJy beam 1 . The individual components fitted to the im- 
ages with JMFIT are summarised in Table [3j where we list the 
flux density, relative position, and size of each component. The 
last two columns show the deconvolved size of the components, 
which describes to a certain degree the intrinsic size of the com- 
ponent, deconvolved from the synthesized beam size. We also 
list the Gaussian components fitted to the uv data (modelfit), 
which are not convolved with any beam. Most of the modelfit 
components are circular instead of elliptical. We plot these com- 
ponents in the last row of Fig.[3]convolved with an artificial cir- 
cular beam with an area equal to the average synthesized beam 
of the rest of the images. As a guide, we plot a line towards the 
direction of the main extended component, which corresponds 
to P.A. of -75, 123, -61, -52, and -74°, respectively. 

All images show a main core component and extended emis- 
sion with changing P.A. During run A, which was obtained soon 



after periastron (orbital phase ~ 0), the source is extended west- 
wards, although the core is better described with two compo- 
nents. This double core has an important contribution eastwards, 
at a distance of ~ 0.5 mas. For the image obtained 24 h later, 
at orbital phase 0.29, a fast morphological change occurs and 
the main extended component is detected towards South-East. 
Run C, obtained soon after apastron, displays bipolar extended 
emission. The emission towards North- West is in fact composed 
by two components, although only one component could be fit- 
ted. The faint component towards South-East could not be fit- 
ted to the data and therefore it is not listed in Table [3] Run D 
is well described with a main core and a bright component to- 
wards North- West. Finally, the images from run E, also obtained 
short after the periastron passage, recover the same morphology 
of run A. This shows that the morphology of LS 5039 is periodic 
on consecutive orbital cycles. 

We also divided each data set in two blocks to check for 
possible intraday variations. We obtained modelfit solutions and 
produced images for the divided uv data sets. However, it was not 
possible to obtain reliable morphological differences because of 
the very different uv coverage at the beginning and the end of 
each run. 

5. Compilation of VLBI observations 

To try to confirm the repeatability of the changing morphologi- 
cal structure of LS 5039 found in our five consecutive days cam- 
paign (project BR127), we compiled the available archival VLBI 
data of the source that provides an angular resolution similar to 
the one explored here. We compiled images from other eleven 
observations conducted at frequencies between 1.7 and 8.5 GHz 
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Table 2. Relative astrometry of the five runs of project BR127. The displacements of the check source J1818— 1 108 are used to 
determine the astrometric uncertainties (see Sect. 14. II ). The mean synthesized beam for J1818— 1 108 in the direct phase referencing 
is 6.5 x 3.3 mas 2 at RA. of 6°, and for the inverse phase referencing it is 6.9 x 2.7 mas 2 at RA. of 7°. We show, for LS 5039 and 
J1825-1718, the beam (HPBW) size and RA. and the relative displacements measured from the average position. These results are 
plotted in Fig. [2] 
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Fig. 3. VLBA images of LS 5039 at 5 GHz from project BR127, obtained during five consecutive days in July 2007. Each column, 
labelled with the run name and orbital phase, corresponds to one epoch (see Table [T]). For each epoch, the first and second rows 
correspond to the self-calibrated images obtained with a natural and uniform weighting scheme, respectively. The third row shows 
the uv components fitted to the data, convolved with a Gaussian circular beam with an area equal to the average synthesized beam of 
all images. The lines show the direction of the main extended component. The restoring beams are plotted in the bottom-left corner 
of each panel. Dashed contours are plotted at -3 times the rms noise of each image and solid contours start at 3 times the rms and 
scale with 2^ 2 . The Gaussian components fitted to the images are listed in Table[3] 



between 1999 and 2009. A log of these observations is shown in 
Table|U T he orbital phases were computed with the ephemerides 
from lCasares et al. N2012al) and have uncertainties below 0.022. 

There are three observations at 5 GHz (6 cm wavelength) 
from 1999 and 2000 that correspond to the projects BP051 and 
GR021A-B, which included the VLBA and the phased VLA ar- 



ray as an a dditional VLB ! statio n. T hese results were already 
presented in iParedes et alJ d2000h and iRibo et a D (120081) . where 
the interested reader can find detailed information on the data 
and the data reduction process. 

Between 2004 and 2006, three runs were observed with the 
VLBA at 8.5 GHz (3.6 cm wavelength) as part of a long-term as- 
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Table 3. Gaussian components fitted to the images shown in Fig. [3] The first two blocks show the fits to the images with natural 
and uniform weighting schemes, respectively, whereas the third block, labelled modelfit, shows fits to the uv data. Each component 
is labelled with the corresponding letter of the epoch (A-E). Columns 2 and 3 list the peak and integrated flux density of each 
component. Columns 4-8 list the position offset and the size fitted to each component. The last two columns list the deconvolved 
size of the components fitted to the images. The modelfit block contains empty values because this fit has no dependence on a 
synthesized beam. 
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Circular 






E2 




9.3 ± 0.5 


-0.00 ± 0.01 


-0.06 ± 0.02 


0.18 ±0.15 


Circular 






E3 




10.5 ± 1.0 


-1.23 ±0.13 


0.36 ± 0.05 


3.66 ±0.16 


Circular 







trometric project (PI: V. Dhawan). The observational codes are 
BD087G, BD105A, and BD105G. The data were correlated with 
a lower sensitivity provided by 128 Mbps, and the final resolu- 
tion is slightly better than for the other projects thanks to the 
higher frequency of the observations. We produced three self- 
calibrated images with natural weight. The obtained synthesized 
beams and noises are show n in Table |H More de tails on the data 
reduction can be found in iMoldon et alJ (120121) . The total flux 
densities recovered for BD087G, BD105A, and BD105G were 
18.6, 17.6, 20.2 mJy, respectively. 

We also present three self-calibrated images from project 
EF018 (PI: R. Fender), which consists of 3 runs (A-C) at 



5.0 GHz obtained with the European VLBI Network (EVN) in 
2007. The runs were separated by 2 days. The antennas used 
were Effelsberg, Westerbork, Jodrell Bank (Lovell), Onsala, 
Medicina, Noto, Torun, Nanshan (Urumqi), Sheshan (Shanghai), 
and Hartebeesthoek. This array provides good sensitivity, al- 
though each run lasted only 4 hours and the elevation of the 
source was lower than in the previous projects. The longest base- 
lines in these observations were from the stations in Nanshan 
(Urumqi) and Hartebeesthoek, although the uv coverage for 
those baselines is poor and the final resolution is lower. The to- 
tal flux densities measured were 26.2, 19.1, and 27.4 mJy, re- 
spectively. The lower flux density and noise of epoch EF018B is 
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Table 4. Log of the VLBI observations of LS 5039. Project 
codes starting with B correspond to VLBA observations, E to 
EVN observations, and G to Global (VLBA+EVN) observa- 
tions. 



EM074 



0.46 



Project 


Epoch 


Orbital 


Freq 


HPBW 




rms 




Y-M-D 


phase 


[GHz] 


[mas - at °] 


[mjy hr 1 ] 


BP051 


1999-05-08 


0.16 


5.0 


3.4 x 1.2 at 


0.4 


0.11 


ER011 


2000-03-01 


0.39 


5.0 


7.6 x 7.0 at 


-14 


0.10 


GR021A 


2000-06-03 


0.49 


5.0 


3.4 x 1.2 at 


0.0 


0.08 


GR021B 


2000-06-08 


0.77 


5.0 


3.4 x 1.2 at 


0.0 


0.11 


BD087G 


2004-06-11 


0.58 


8.5 


2.6 x 1.1 at 12.2 


0.11 


BD105A 


2005-06-11 


0.02 


8.4 


2.9 x 1.3 at 


7.9 


0.09 


BD105G 


2006-01-29 


0.51 


8.4 


3.4 x 1.2 at 


7.3 


0.10 


EF018A 


2007-03-01 


0.89 


5.0 


5.9x4.2 at 


5.3 


0.14 


EF018B 


2007-03-03 


0.40 


5.0 


6.9x4.1 at 


7.3 


0.04 


EF018C 


2007-03-05 


0.91 


5.0 


4.5 x 3.6 at 


8.1 


0.11 


BR127A 


2007-07-05 


0.03 


5.0 


6.0 x 2.4 at 


-0.4 


0.06 


BR127B 


2007-07-06 


0.29 


5.0 


6.4 x 2.4 at 


-3.0 


0.06 


BR127C 


2007-07-07 


0.55 


5.0 


6.0 x 2.3 at 


-2.2 


0.06 


BR127D 


2007-07-08 


0.80 


5.0 


5.9 x 2.2 at 


-2.3 


0.07 


BR127E 


2007-07-09 


0.06 


5.0 


6.3 x 2.4 at 


-3.2 


0.07 


EM074 


2009-03-07 


0.46 


1.7 


26 x 5.8 at 


7.3 


0.04 




20 10 -10 -20 -30 -40 
R.A. offset [mas] 



Fig. 5. Self-calibrated image of LS 5039 at 1.7 GHz obtained 
with the EVN in 2009. The plot parameters are the same as in 
Fig- HI although a larger angular scale is shown here. 



probably an artifact of the amplitude self-calibration and not a 
real change in the source flux density. 

In Fig. 3] we show the self-calibrated images of these 
projects. The images in the panels are ordered in increasing 
phase, from left to right and from top to bottom. We distributed 
them in four groups, which correspond to images with similar or- 
bital phases and that show approximately a similar morphology. 
For a quick reference, the inset table included in Fig. [4] shows 
the dates of the images in the panels. Although they were ob- 
tained with different instruments, frequencies, and at very differ- 
ent epochs, the images in each group share common morpholog- 
ical trends. Each group can be approximately characterised with 
the morphology described in Sect. l4.2l for the runs BR127A-D. 

5.1. Morphology at larger scales 

We observed LS 5039 at 1.7 GHz (18 cm wavelength) with 
the EVN in 2009, project code EM074. This is a deep pointed 
observation on LS 5039, conducted without phase referenc- 
ing. The total time on source was 5.8 h, and 9 stations par- 
ticipated: Effelsberg, Westerbork, Cambridge (32 m), Jodrell 
Bank (Lovell), Onsala, Medicina, Noto, Torun, and Urumqi. At 
this frequency the resolution is lower, although the long East- 
West baselines provide good resolution in right ascension. We 
produced a uniform weighting scheme image with a synthe- 
sized beam of 26 x 5.8 mas 2 at P.A. 7° and a rms noise level 
of 0.04 mjy beam -1 . The total flux density of the source was 
33.5 mjy at 1.7 GHz. The self-calibrated image is shown in 
Fig. |5] The source shows a morphology more similar to the im- 
ages in group 2, although it was observed at orbital phase 0.46, 
which is closer to the ones of group 3. This can be explained by 
the longer lifetime of electrons emitting synchrotron radiation at 
lower frequencies. 

On the other hand, the self-calibrated image from project 
ER011 provides lower angular resolution than the images in 



Fig. |H a nd it is not reproduc ed here because it is already pub- 
lished in Pared eTet al.l J2002). The image was obtained at orbital 
phase 0.39. The innermost region of the core shows extended 
emission towards South-East, as it happens in the group 2 images 
in Fig.|4] although additional bipolar diffuse components ex tend 
up to ~ 20 mas in opposite directions. iParedes et al.l J2002) also 
present a lower resolution image obtained with MERLIN show- 
ing extended emission with P.A. of 150° up to 170 mas. Finally, 
we note that lower resolution VLA images show some hints of 
extended emission in the same North- West/South -East direction 
at angular distances of > 0.1" (M artf et all 1998). This confirms 
that LS 5039 contains an additional extended and diffuse com- 
ponent that cannot be traced with the high-resolution images de- 
scribed in Sect. l4.2l andl5l 

6. Model 

In this section we try to reproduce the changing radio morphol- 
ogy of LS 5039 by modelling the emission produced by an out- 
flow of electrons accelerated as a consequence of the interaction 
of the stellar wind and the wind of the putative young pulsar. 
We are interested in describing the orientation of the P.A. of the 
emission up to scales of ~ 5-10 mas. The formation and evo- 
lution of the outflow is complex, and hydrodynamical and mag- 
netohydrodynamical simulations are needed to describe in detail 
the source radio emission (see the discussion in Sect.[TJi. Our aim 
here is to check if the general scenario of colliding winds is com- 
patible with the main features and the continuous morphological 
changes found at scales up to 10 mas. 

We assume a mass for the compa nion star of 33 M H a nd an 
eccentricity of the system of 0.35 dCasares et alj|2012a ). The 
semimajor axis of the orbit for a 1.4 M Q neutron star is ap- 
proximately a = 2.4 x 10 12 cm, or 0.05 ma s for a distance 
to the source of 2.9 kpc dCasares et al.ll2012ah . The semimajor 
axis does not change significantly for compact object masses be- 
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BD105A 




R.A. offset [mas 



Project 


Obs. Date 


BD105A 


2005-06-11 


BR127A 


2007-07-05 


BR127E 


2007-07-09 


BP051 


1999-05-08 


BRI27B 


2007-07-06 


EF018B 


2007-03-03 


GR021A 


2000-06-03 


BD105G 


2006-01-29 


BRI27C 


2007-07-07 


BD087G 


2004-06-11 


GR02IB 


2000-06-08 


BR127D 


2007-07-08 


EF018A 


2007-03-01 


EF018C 


2007-03-05 



BD087G 



0orb = 0-58 




EF018C 



: 0.91 



1 1 1 




i i i 


5.0 GHz 
i i 



10 



5 0-5 
R.A. offset [mas] 



-10 



Fig. 4. VLBI self-calibrated images of LS 5039 at 5.0 and 8.5 GHz (see Table|4]for details on the projects). The plotting parameters 
are the same as in Fig. [3] The panel order is in increasing orbital phase from left to right and from top to bottom, and the images 
have been grouped according to similar morphological features. The project code and orbital phase are quoted on the top part of 
each panel. The images from project BR127, obtained during the same orbital cycle, have bold axes. Dashed contours are plotted 
at -3 times the rms noise of each image and solid contours start at 3 times the rms and scale with 2^ 2 , except for GR021A-B and 
BP051, which start at 5cr as in the original publications. The rms noise of each image can be found in Table|4] 



low 5 M Q . iBosch-Ramon et al" identified regions where 

strong shocks are produced at distances between 5 and 10a from 
the massive star, towards the direction of the pulsar. Here we will 
assume that the extended radio emission is produced by elec- 
trons accelerated at a distance of 10 times the semimajor axis 



of the orbit, ~0.5 mas at 2.9 kpc, and neglect the high-energy 
processes taking place at shorter distances. We show a sketch 
of the scenario in the first row of panels of Fig. |6j where the 
grey area indicates the flow of particles accelerated at a dis- 
tance of 10a from the star and are travelling away from it. We 



10 



Moldon et al.: Periodic morphological changes in the radio structure of LS 5039 



consider that the injected electrons follow an energy power law- 
distribution, N e oc yj p . The glob al radio spectral index of the 
source is -0.46 (M arti et alj!9 98). and therefore we use an elec- 
tron index p — 2. For simplicity, we assume a constant injec- 
tion rate of 2inj = 10 35 erg s _I (see below) along the orbit, for 
electron energies of 1 < y e < 10 4 . Higher electron energies 
do not significantly contribute to the radio emission at the con- 
sidered frequencies. On the other hand, a range for the lower 
electron energy of ~ 1-100 slightly modifies the most inner 
part of the emission, but does not significantly affect the larger 
source structure at 2-5 mas. The model is not detailed enough to 
constrain this value and therefore we assume the minimum en- 
ergy possible. We note that much higher minimum energies are 
expected in the primary shock bet ween the stellar and the pul- 
sar wind dKennel & Coronitil [l984) but we are here consid ering 
a secondary shock at larger distances. Zab alza et al.l d201 lb in- 
ferred a spin-down luminosity for the putative pulsar in LS 5039 
of L s d ~ (3-6) x 10 36 erg s _1 , and therefore our injected energy 
accounts for ~ 2% of this pulsar spin-down luminosity. The en- 
ergy spectrum was divided in 200 logarithmic intervals for the 
computations. 

For distances from the massive star larger than 10a, the evo- 
lution of the electron distribution is determined by the radia- 
tive losses (synchrotron and inverse Compton in the Thomson 
regime) and adiabatic cooling. We assume that the particles 
travel away from the massive star forming an outflow with con- 
ical shape (see Fig. [6j. This is conical shape is a coarse esti- 
mation, as the ambient pressure and the hy drodynamical insta- 
bilitie s severely modify the flow structure (Bosch-Ramon et all 
120121) . We use an opening angle for the cone of 20°, of 
the same order as the post-shoc k flow structure seen in 
iBosch-Ramon et all (l2012h (see also lBogovalov et al ] l2008h . al- 
though the final flux density distribution is not very sensitive 
to this value. We note that this angle is much lower than the 
asymptotic opening angle of the contact discontinuity, which for 
this system is about ~ 75° assuming a pulsar with spin-down 
luminosity of 5 x 10 36 erg s _1 , with a stellar wind with a mass 
loss rate of 2.5 x 10~ 7 M H yr~' and w ind velocity at infinity of 
2.4x 10 3 km s" 1 dMcSwain et al.l201 lh . We compute the electron 
energy distribution density n(r, y), where r is the linear distance 
from the massive star, taking into account the conservation of 
the number of particles along the outflow through the continuity 
equation N{r, y)dy = N(ro, yo)dyo. 

The outflow is bent as the particles travel away from the sys- 
tem because of the orbital motion of the pulsar (see second row 
of panels of Fig.|6]l. The emitting region is bent in the opposite 
directi on of the orb ital motion of the pulsar and it forms a spi- 
ral (see lDubuil2006l for a general description). The flow, similar 
to a cometary tail, is a projection of the orbit assuming a bal- 
listic motion. However, the outflow, and in general the whole 
interaction structure, in fact should not follow a ballistic mo- 
tion because of the effect of the Coriolis force and consequently 
the resultin g spiral should be less open than for simple ballis- 
tic m otion (Bosch-Ramon & Barkovl l201 U IBosch-Ramon et al. 
120121) . Nevertheless, we approximate the outflow shape assum- 
ing a constant advection velocity of 0.04c, which traces an struc- 
ture similar but more o pen than the one found in the hydrody- 
namical simulations bv lBosch-Ramon et alJ (1201 2l) . because we 
find a better match with the images in this case. We note that 
those simulations were computed for a slightly different system 
and using a circular orbit. 

We compute the electron density along the flow for r be- 
tween 10a and 140a, although the farthest part (> 100a) does 
not contribute significantly to the emission. We divide the flow 



in 200 steps in r. The sizes of the steps is determined by the 
displacement of the pulsar along the orbit and are computed at 
constant increments of the true anomaly of the pulsar. Therefore, 
steps close to periastron are shorter, corresponding to time in- 
tervals of about 500 s each, whereas steps close to apastron 
are longer, corresponding to time intervals of about 3000 s 
each. These intervals are always at least one order of magni- 
tude shorter than the fastest cooling time at each r and y e . For 
the lower energy electrons, adiabatic cooling dominates along 
the flow, whereas for more energetic electrons, inverse Compton 
losses dominate up to distances of < 4 mas. We compute the 
emissivity of each section of the flow assuming a magnetic field 
of B = 1.0 [ 2 - 4 x'o' 2cm ] g, which corresponds to B ~ 0.1 G in the 
acceleration region. We note that this magnetic field is similar 
to the equipartition magnetic field inferred for no n-thermal par- 
ticles from VLBI observations, which is ~ 0.2 G dParedes et al.l 
l2000l) . 

The observed outflow geometry strongly depends on the ori- 
entation of the orbit on the sky, which is basically determined 
by two angles: the longitude of the ascending node, Q, mea- 
sured from North to the East, and the inclination of the orbit, 
2, which is a rotation with respect to the node axis (the orbital 
elements are defined in lSmartl[l930l) . In Fig. [6] we show different 
projections of these two angles at orbital phase 0.5. The images 
in Fig. [4] show a privileged North-East/South- West direction, in 
a RA. of ~ 130°. This constant direction for most part of the 
orbital phas es suggests a high inclination of the orbit, already 
discussed in iRibo et al.1 (12008b . The orbital elements determine 
the flow orientation and the final flux density distribution. We 
searched for the combination of Q and i that better describes the 
source structure seen in Fig. [4] We used the electron density and 
the system geometry to compute the synchrotron emissivity of 
the electrons and project it on the sky. For every orbital phase 
in Fig. H] we computed the sky flux density distribution at the 
corresponding frequency and convolved it with a beam equal to 
the corresponding synthesized beam of the image. The model 
accounts for part of the core component and the extended emis- 
sion. We produced synthetic images in steps of 5° in i and Q. 
independently and we found that the images in Fig. [4] are bet- 
ter reproduced by the model for an inclination of i ~ 70° and 
Q ~ 130°. Based on visual inspection, a reasonable range for 
these parameters is i = 60 - 75° and Q = 120 - 140°. The 
rotation of the pulsar a round the sta r is counter-clockwise, and 
therefore i < 90° (see ISmartl [19301) . Soon after the periastron 
passage the pulsar transits behind the star plane and after apas- 
tron it is in front of the star plane. For the inclination obtained 
the pulsar position is eclipsed during its superior conjunction. 

The resulting synthetic images and the projected orbit are 
shown in Fig. [7] The first contour level for all plots is set to 
0.06 mJy beam 1 . The model produces total flux densities of 
around 10-16 mJy, which accounts for the extended flux den- 
sity recovered with VLBI. However, this model does not account 
for the emission produced at distances below < 0.5 mas, so we 
artificially included an additional 10 mJy point-like source at 
the position of the binary system to visually help to compare 
the results with the images in Fig. |4] This component approxi- 
mately accounts for the peak flux density of the core component 
of the images. With this ad-hoc component, we recover a total 
flux density of the order of 15-18 mJy at 5 GHz and 13-16 mJy 
at 8.5 GHz, for the considered orbital phases, which are slightly 
below the measured values. 

This model describes correctly the main features in most of 
the observed images, except for the observations shortly after 
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the periastron passage. For that case (group 1), our model only 
accounts for the secondary component eastwards, but not for the 
main extended component towards West. The rest of the images 
are well described by this simple model. It produces a bright 
extended component for orbital phases 0.15-0.40 (group 2) and 
continuously develops extended emission towards North- West 
between phases 0.5 and 0.6 (group 3), while preserving the bipo- 
lar structure. The extended emission is then dominated by the 
North-West component until the next periastron passage. We 
note that this model predicts bipolar extended emission at cer- 
tain orbital phases, although the outflow is not bipolar itself. In 
particular, the projected flow at phase 0.5, shown in the bottom- 
right panel of Fig. [6] accounts for the bipolar flux density dis- 
tribution seen for instance in the image of project GR021A (see 
Figs. |4] and |7). At the same time, this approach justifies that the 
bipolar emission is not complet ely symmetric an d that the two 
components show different P.A. (Ribo et al. 2008]). 

The differences between the model (Fig. |7]l and the real 
data (Fig.|4]l, in particular the missing component in the images 
from group 1, can be explained by the oversimplification of the 
flow structure and the constant physical conditions at all orbital 
phases. For instance, we assumed that the acceleration region 
was at a constant distance from the massive star and that the en- 
ergy injection was constant, although it is expected that these 
parameters depend on the orbital phase in an eccentric binary 
system. Therefore, additional electron populations, as well as 
more realistic time-dependent physical conditions should be in- 
cluded to produce this additional component. The model should 
also include considerations on the relativistic Doppl er boostin g 
possibly affecting the synchrotron emission dDubus et al.l l2010). 
However, including additional particle populations and assum- 
ing time-dependent physical conditions would require the addi- 
tion of several free parameters, which is outside of the scope of 
this paper. Therefore, we do not try to include a complete flux 
density analysis or to predict flux density distributions. Also, 
predictions of displacements of the peak of the emission, such as 
the hint found in Sect. 14.11 cannot be obtained by this model, as 
these are produced at smaller distances. Despite this caveats, this 
model traces, except for the phases close to the periastron pas- 
sage, all the main features of the VLBI images and the morpho- 
logical variability measured in multifrequency and multiepoch 
observations, even when assuming simple and constant physical 
conditions. 

7. Discussion and conclusions 

VLBA observations at 5 GHz during five consecutive days show 
that LS 5039 displays orbital morphological variability, showing 
one sided and bipolar structures, but recovering the same mor- 
phology when observing at the same orbital phase. The P.A. of 
the extended emission with respect to the core changes signifi- 
cantly on timescales of one day. Although the total flux density 
remains approximately constant during the orbital cycle, the flux 
density of the different components varies. Remarkably, the RA. 
of the extended emission changes by ~ 160° between phases 
0.03 and 0.29 (24 h). The sensitivity of the observations is not 
enough to study changes on shorter timescales, within a few 
hours. We measured a displacement of the peak of the emis- 
sion between orbital phases 0.29 and 0.55 of 2.3 + 0.6 mas in 
RA. 154°, nearly in the same direction of the extended emis- 
sion. However, we consider this measurement as a hint of peak 
displacement because, given the limitations of the astrometry 
of these observations, we would requ ire a 5- c r sign ificance to 
consider it a significant displacement. IDubusI (|2006) predicted 



a continuous peak displacement forming an ellipse of a few 
mas along the orbit. The peak of the emission of the gamma- 
ray binary LS I +61 303 at 8.4 GHz traces an ellipse of semi- 
major axis 1.5 AU, about four times the binary semimajor axis 
dDhawan et alj|200"6h . Unfortunately our astrometric accuracy is 
not enough to measure these displacements. The source mor- 
phology during two consecutive periastron passages (run A and 
E) is very similar, showing that the changes within the same or- 
bital cycle are periodic. An interesting future project would be 
aimed to obtain accurate astrometry to explore the displacements 
of the peak of the emission, in particular at different frequen- 
cies, and also to measure variations of the source structure on 
timescales of a few hours. 

We analysed the available VLBI observations of LS 5039 
in a consistent way. Images from data obtained between 1999 
and 2009 at frequencies between 1.7 and 8.5 GHz show that 
the morphology at similar orbital phases is similar. When or- 
dered in increasing phase, the images show an approximately 
continuous change. The observational conclusion of our VLBI 
analysis is that the morphological changes of LS 5039 are pe- 
riodic and show orbital modulation stable over several years. 
The gam ma-ray binary LS 1 +61 303 also shows a similar be- 
haviour (Dha wan et alj 12006). Therefore, all gamma-ray bina- 
ries are expected to display periodic orbital modulation of their 
VLBI structure. 

A simple model of an expanding outflow of relativistic elec- 
trons accelerated at a distance of ~ 10a and travelling away from 
the system can account for the main changes seen in the extended 
emission of LS 5039. The changes are explained for an inclina- 
tion of the orbit i ~ 70° and a longitude of the ascending node 
of Q. ~ 130°. In the last row of panels in Fig.|6]we show the ori- 
entation of the orbit on the sky projected with these parameters. 
The inclination of the orbit has deep implications in the physi- 
cal properties of the binary system, because the orbital param- 
eters restrict the mass of the compact object within the system. 
In particular, a ssuming a mass fun ction of the binary system of 
0.0045(6) M dCasaresetal.ll2012al) . a mass of the star of 33 M Q 
implies that the mass of the compact object is 1 .8-2.0 M Q , for an 
inclination between 75 and 60°. However, the mass of the main 
star is barely constrained between 20 and 50 M , and conse- 
quently the inclinations inferred from the model are compatible 
with masses between 1.3-1.5 and 2.4-2.7 M , for the lower and 
the higher stellar mass limits, respectively. The uncertainty in 
the mass function adds an additional uncertainty of ~ 0.1 M . 
The values obtained are in any case below 3 M , and therefore 
are compatible with the presence of a neutron star in the system. 
However, we note that this was one hypothesis of the model and 
therefore cannot be a proof. After applying a simple model, we 
conclude that that the presence of a young non-accreting pulsar 
in LS 5039 agrees with the observed source structure at different 
orbital phases and produces a coherent scenario. However, this 
model does not account for all the emission produced at shorter 
distances and also fails to explain the main component of the ex- 
tended emission seen shortly after the periastron passage. This 
is possibly due to unaccounted variability in the injection rate, 
geometry and localisation of the acceleration zone, and flow ve- 
locity of the electrons. 

For gamma-ray binaries, the nature of the compact object, 
either accreting or not, is a fundamental ingredient for any 
model aimed to understand these peculiar systems. The most 
straightforward indication would be the detection of pulsations. 
However, direct detection may be unfeasible in binary systems 
with close orbits and powerful ma ssive companions because of 
free-free absorption (Dubus 2006). In addition, the mass func- 
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tion of the known gamma-ray binaries does not allow to dis- 
criminate between black holes or neutron stars. Therefore, an 
ongoing debate is still open regarding the nature of the compact 
object of gamma-ray binaries. We have presented high resolution 
radio images of LS 5039 with which the scenarios can be tested. 
A first approach has shown that the source radio morphology can 
be explained with the presence of a young non-accreting pulsar. 
PSR B 1259-63, the only gamma-ray binary wit h a confirmed 
pulsar, shows variable extended radio structures (Moldon et al. 
1201 lah . and LS I +61 303 shows var iable and periodic mor- 
phological changes dDhawan et alj|200"6h . The common features 
among these three systems suggest that the known gamma-ray 
binaries contain young non-accreting pulsars. 
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Fig. 6. Sketch of the orbit orientation and the emission region of the model. From left to right the panels show three different scales: 
the orbit (a), the acceleration region (10a), and the whole flow (100a). The particle acceleration occurs at a distance of 10a from 
the massive star (blue circle) behind the pulsar (empty green circle) and the flow of emitting particles (grey area) expands forming 
a cone while travelling away from the binary system. The star plane is the plane of the sky intersecting the star. The colours in the 
orbit and the flow indicate regions situated in front of the star plane (red, closer to the observer than the star) and behind the plane 
(yellow, farther than the star). The curved arrows indicate the orbital motion of the pulsar. The first row shows the general scenario 
without orbital motion, whereas the rest of the rows include the orbital motion and different projections of the orbit, indicated in the 
labels. The last row of panels shows the orbit projected with our best estimate of i and £1 
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p mb = 0.02 <f, mb = 0.03 <f> olb = 0.06 
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Fig. 7. Synthetic images produced with our model, computed for the corresponding orbital phases of the different VLBI projects 
and convolved with the corresponding beam shown in Fig. [4] The black line traces the projected position of the flow axis of cooling 
particles, which is computed for linear distances from the star between 1.6 and 22 AU. The crosses mark the position of the main 
star, set at the origin. The panel in the top-right corner shows the system orbit, projected with an inclination of 70° and a longitude 
of the ascending node of 130°. 
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